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ABSTRACT 



We investigate differential rotation in rapidly rotating solar-type stars by 
means of an axisymmetric mean field model that was previously applied to the 
sun. This allows us to calculate the latitudinal entropy gradient with a rea- 
sonable physical basis. Our conclusions are as follows: (1) Differential rotation 
approaches the Taylor-Proudman state when stellar rotation is faster than so- 
lar rotation. (2) Entropy gradient generated by the attached subadiabatic layer 
beneath the convection zone becomes relatively small with a large stellar angu- 
lar velocity. (3) Turbulent viscosity and turbulent angular momentum transport 
determine the spatial difference of angular velocity AQ. (4) The results of our 
mean field model can explain observations of stellar differential rotation. 

Subject headings: Sun: interior — Sun: rotation — Stars: interior 
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INTRODUCTION 



Our sun has an eleven-year magnetic activity cycle, which is thought to be sustained 
by the dynamo mot ion of interna l ionized plasma, i.e., a transformation of kinetic energy 



to magnetic energy (IParker 



19551 ). Our understanding of the solar dynamo has significantly 



improved during the past fifty years, and some kinematic studies can now reproduce solar 



magnetic features suc h as equatorward migration of sunspots and pq 



of the magnetic field (iChqudhuri et al. 



2001 



Charbonnea 



. 



2005 



1995 



Hotta fc Yokoyama 



3ik pati fc Charbonneau 



2010a 



eward migration 



1999; 



Kuker et al. 



bj). The most important mechanism of 



the solar dynamo is the Q effect, the bending of pre-existing poloidal magnetic fields by 
differential rotation and the generation of toroidal magnetic fields. Thus, the distribution of 
the differential rotation in the convection zone is a significant factor for the solar dynamo. 
Using helioseismology, it has recently been shown t hat the solar internal d ifferential rotation 



Thompson et al. 



20031 ). meaning the 



is in a non- Taylor- Proudman state (see review by 
iso-rotation surfaces are not parallel to the axis. 

Based on solar observations, it is known that Ca H-K fluxes can be a signature of stellar 
chromos pheric a c tivity , and such chromospheric signatu res are in correlation with magnetic 



activity. 



Wilson! ( 11968 



1978|) and 



Baliunas et al. 



( 119951 ) discuss a class of stars that shows a 



periodic variation in Ca H-K fluxes, which suggests that they have a magnetic cycle similar 
to our sun. It is natural to conjecture that such magnetic activity is maintained by dynamo 
action. Various studies have been conducted to investigate the relationship between stellar 



angular velocity f2 and its latitudinal difference AQ i.e., AQ oc fffi, w 



range of n is < n < 1 (IDonahue et al. 



1996 



Reiners fc Schmitt 



2003 



rere the suggestec 



Barnes et al. 



20051). 



This means that the angular velocity difference AQ increases and the relative difference 
AQ/Q decreases with increases in the stellar rotation rate Qq. 



In this paper, we investigate differential rotation in rapidly rotating stars using a mean 
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field framework. Our study is based on the work of iRempell (J2005bl ). in which he suggests 
the importance of the role of the subadiabatic layer below the convection zone in order to 
maintain a non-Taylor-Proudman state in the Sun. The aim of this paper is to use a mean 
field model to analyze firstly the dependence of the morphology of differential rotation on 
stellar angular velocity, and secondly the physical process which determines the observable 
angular velocity difference AQ. According to our knowledge, this is the first work which 
systematically discusses the application of Rempel's (2005b) solar model to stars. 



Other research adopts another approach to the use of mean field models 



1995 



or the analysis 



Kuker fc Stix 



of dif ferential rotation in rapidly-rotating stars (jKitchatinov fc Riidiger 
20011 ). In these studies, the non-Taylor-Proudman state is sustained by anisotropy of 
turbulent thermal conduction. This anisotropy is generated by the effects of stellar rotation 
on convective turbulence. 

Three-dimensional numerical studies on stellar differential rotation also exist 



( iBrown et al. 



2008; 



Miesch fc Toomre 



20091 ) . In these studies, they resolve stellar thermal 
driven convection and can calculate a self-consistent turbulent angular momentum transport 
and anisotropy of turbulent thermal conductivity. The subadiabatic layer below the 
convection zone, however, is not included. The effects of anisotropy of turbulent thermal 
conductivity and the subadiabatic layer are discussed in this paper. 



2. MODEL 

Using numerical settings similar to those of Rempel's (2005b), we solve the axisymmetric 
hydrodynamic equations in spherical geometry (r,6), where r is the radius, and 9 is the 
colatitude. The basic assumptions are as follows. 



1. A mean field approximation is adopted. All processes on the convective scale 
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are parameterized. Thus, the coefficients for turbulent viscosity, turbulent heat 
conductivity, and turbulent angular momentum transport are explicitly given in the 
equations. 

2. The perturbations of the density and pressure associated with differential rotation are 
small, i.e., p\ <C po an d p\ <C po- Here p and po denote the reference state density 
and pressure respectively, whereas p\ and p\ are the perturbations. We neglect the 
second-order terms of these quantities. Note that the perturbation of angular velocity 
(fii) and meridional flow (v r , Vg) are not small. 

3. Since the reference state is assumed to be in an energy flux balance, the entropy 
equation includes only perturbations. 

2.1. Equations 

We do not use the anelastic approximation here. The equations in an inertial frame 
can be expressed as 

dpi 1 <9 2 1 d . 

at H or r sin 9 o9 



dv r 8v r Vg 8v r Vq 1 

dt dr r d9 r p 



dpi 
Pl9+ ^r~ 



+ (2^^ + nl)r sin 2 9 + — , (2) 

Po 



dv e dv e Vg dvg V r V 9 1 1 d^ f 2 Fg 

— = -v r — — — - + (2fi fi 1 + Q 1 )rsm#cos# + — , (3) 

at or r at) r p$r 06 po 

^r = --2^ r2 ^ + n 01 ~ ^r-nU^ 2 e ^ + ^ + — Fj H> ( 4 ) 

at r ^Q r l rsm 39 p r smv 

dsi dsi vgdsi 7<5 7-1 1 

7T = -^ - 7^ + ^ + ^T g + ^ dlv (^oT grad Sl ), (5) 

where Qq is a constant value that represents the angular velocity of the rigidly rotating 
radiative zone. We set it as a parameter in Table [TJ 7 is the ratio of specific heats, with the 
value for an ideal gas being 7 = 5/3. K t is the coefficient of turbulent thermal conductivity. 
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5 = V — V a d represents superadiabaticity, where V = d(\nT)/d(\np) (see §2.21) . g denotes 
gravitational acceleration. Following from this, the perturbation of pressure p\ and pressure 
scale height H p are expressed as 



. P 1 
Pi=Po\ 7 1- si 

Po 
Po 



H n 



Po9 



(6) 
(7) 



Si is dimensionless entropy normalized by the specific heat capacity at constant volume c v . 
Turbulent viscous force F follows from 



}_d_, 2r _L 

r 2 dr r sin 9 09 



9 > ■ a d > Ree + R<t 



r 2 R r e) + 



1 9 
with the Reynolds stress tensor 



r sin 9 89 



d , 

'sin QRee) + 







r sin <9# 



sin 





r 


i 


R r o 


— R<f>4> 


cotfl 




r 




lLptp 


+ Re<f> 


cot# 



+ 



(8) 



(9) 



(10) 



R, 



ik — PO 



vu 



I E ik - -<5 ifc divv j 



^tiA 



ik 



;n) 



iere Vt v is the coefficient of t urbulent viscosity and u t \ is the coefficient of the A effect 



Kitchatinov &: Riidiger 



19951 ). a non-diffusive angular momentum transport caused by 
turbulence. u tv and z/ tl are expected to have the same value, since both effects are caused 
by turbulence, i.e., thermal driven convection. We discuss this in more detail in §2.31 En, 
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denotes the deformation tensor, which is given in spherical coordinates by 

E =2^ 
or 



E r g 

E r( f) 

Eqa, 



1 dv g v r 

r ou r 

2. 

- [V r + Vg COt l 

r 

d_ (vf 

dr V r 



Eg r = r— i — 

E^ r = r sin 6? 
E M = sin 61 



1 dv r 

r 80 



dr 
9Qi 



90 



(12) 
(13) 
(14) 
(15) 
(16) 
(17) 



An expression for the A effect (Ajj.) is given later. The amount of energy that is converted 
by the Reynolds stress from kinematic energy to internal energy is given by 



Q — 2_^ ^EikRik- 



'11 



i.k 



2.2. Background Stratification 



We use an adiabatic hydrostatic stratification for the spherically symmetric reference 
state of Poj Po an d T . Gravitational acceleration is assumed to have ~ r~ 2 dependence, 
since the radiative zone (r < 0.65i? o ) has most of the solar mass. This is expressed as, 



Po(r) = Pbc 



p {r) 
To(r) 
g(r) 



= Pbc 

= ?bc 

fl'bc 



1 + 
1 + 



7 - 1 r b 



i / ^"bc -. 



7 H hc V r 

7 - 1 n> 



< / '"be -. 



7 H hc V r 



1/(7-1) 



7/(7-1) 



7 H hc V r 



x 7 ~ 1 r hc /rv _ ! 



r bc 



(19) 

(20) 
(21) 
(22) 



where pbc, Pbc, ^bc, Hhc = Pbc/(pbcfi'bc) and ^bc denote the values at the base of the 
convection zone r = r^, c of density, pressure, temperature, pressure scale height and 



gravitational acceleration, respectively. In this study we use r^ c = 0.71R Q , with R & 
representing the solar radius (R Q = 7 x 10 10 cm). We adopt solar values pbc = 0.2 g cm" 3 , 
Phc = 6 x 10 13 dyn cm -2 , T bc = mp hc /(k B p hc ) ~ 1.82 x 10 6 K and g hc = 5.2 x 10 4 cm s~ 2 , 
where k-Q is the Boltzmann constant, and m is the mean particle mass. Fig. [1] shows the 
profiles of background density, pressure and temperature, and gravitational acceleration. 

Although the real sun's stratification is not adiabatic in the convection zone, our 
reference state is valid, since the absolute value of superadiabaticity is small. In order to 
include the deviation from adiabatic stratification, we assume superadiabaticity 5 has the 
following profile: 



"conv T „ \0 OS O conv ) 



i I ' ' tran 

tanh 



d 



tran 



(23) 



Here <5 os and 5 conv denote the values of superadiabaticity in the overshoot region. r tran and 
c^tran denote the position and the steepness of the transition toward the subadiabatically 
stratified overshoot region, respectively. Superadiabaticity in convection zone is define as 

^conv = <5c ^ , (24) 

'max 'sub 

where r max denotes the location of the upper boundary. We specify 5 os = —1.5 x 10~ 5 , 
'"tran = O.725_R , r su b = O.8i? and rf tran = <i su b = 0.0125i? Q in our simulations. 5 C is took as 
a free parameter. The entropy gradient can be expressed as 

ds _ jS 

dr~ H p [Zb) 

The third term of eq. (JSJ), v r ^5/H pi includes the effect of deviations from adiabatic 
stratification. The term indicates that an upflow (downflow) can make negative (positive) 
entropy perturbations in the subadiabatically stratified layers (5 < 0). 
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2.3. Diffusivity Profile 



We assume the coefficients of turbulent viscosity and thermal conductivity to be 
constant within the convection zone, and these smoothly connect with the values of the 
overshoot region. We assume that the diffusivities only depend on the radial coordinate: 

r ~ r tran + A 



with 



Vt\ 






^0v 

2 



1 + tanh 



(L 



Mr), 



1 + tanh 



r ~ r tran + A 



«t 



, «0 



1 + tanh 



r ~ r tran + A 



/c(0, 



/e(r), 



fc(r) = - 



1 + tanh 
-l 



A = d KV tanh (2a 



r -r bc 

-1), 



(26) 

(27) 
(28) 



(29) 
(30) 



where z/ 0v , i^oi an d K o are the values of the turbulent diffusivities within the convection zone, 
and u os and k os are the values in the overshoot region. We specify z/ i = Ko\ = 3 x 10 12 cm 2 s _1 , 
^os = 6 x 10 10 cm 2 s _1 and n os = 6 x 10 9 cm 2 s _1 , and we treat u 0v as a parameter. a Kl/ 
specifies the values of the turbulent diffusivities at r = r tran , i.e., v tv = u os + a KU u 0v , 
vt\ = ol kv vq\ and n t = ^ os + ct KU n at r = r tran . d\, c and d RV are the widths of transition. 
We specify a KV = 0.1, db c = 0.0125i? o and d KV = 0.025i? Q . As already mentioned, the 
coefficients for turbulent viscosity and the A effect are different in our model from those 
of Rempel's (2005b). There are two reasons for this. One is that we intend to investigate 
the influence of both effects on stellar differential rotation separately (see §4.2j) . The other 
reason is that the formation of a tachocline in a reasonable amount of time requires a finite 
value (though small) for the coefficient of tu rbulent viscosit y even in the radiative zone, in 



which there is likely to be weak turbulence (iRempel 
z/ tv , z/ tl and K t . 



2005bl ). Fig. |2] shows the profiles of 
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2.4. The A Effect 

In this study we adopt the non-diffusive part of the Reynolds stress, called the A effect. 
The A effect transports angular momentum and generates differential rotation. The A effect 
tensors are expressed as 

A r ^ = A^. = +L(r,9)cos(9 + \), (31) 

A^ = A^ = -L(r,0)sin(0 + A), (32) 

where L(r, 9) is the amplitude of the A effect and A is the inclination of the flux vector with 
respect to the rotational axis. We use for the amplitude of the A effect the expressions 

f(r, 9) = sin' 9 cos 9 tanh [ **,' J , (33) 

L(r,9)=A n f{ ^\ M , (34) 

max|/(r, 9)\ 

where d = 0.025i? Q . A and A are free-parameters. The value of I needs to be equal to or 

larger than 2 to ensure regularity near the pole, so we set I — 2. The A effect does not 

depend on v r , Vg or Qi, meaning it is a stationary effect. We emphasize that the A effect 

depends on stellar angular velocity O , since the A effect is generated by turbulence and 

Coriolis force. The more rapidly the star rotates, the more angular momentum the A effect 

can transport. The dependence of A and A on stellar angular velocity is discussed in §4.31 

2.5. Numerical Settings 



19841 ) . 



Using the modified Lax-Wendroff scheme with TVD artificial viscosity (ID avid 
we solve Equations (JIJ-flHJ) numerically for the northern hemisphere of the meridional plane 
in 0.65-R Q < r < O.93i? and < 9 < %/2. We use a uniform resolution of 200 points in the 
radial direction and 400 points in the latitudinal direction in all of our simulations. Each 
simulation run is conducted until it reaches a stationary state. All the variables pi, v r , vg, 
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Qi and s± are equal to zero in the initial condition. At the top boundary (r = O.93_R ) we 
adopt stress-free boundary conditions for v r , vg and Qi and set the derivative of Si to zero: 



dv r 

dr 
d_ fv_e_ 

dr V r 

dr 
ds\ 

dr 



0. 



0. 
0. 



(35) 
(36) 
(37) 
(38) 



The boundary conditions for v r , vg and si at the lower boundary (r = 0.65-R Q ) are the same 
as those at the top boundary. Differential rotation connects with the rigidly rotating core 
at the lower boundary, so we adopt Qi = there. At both radial boundaries, we set p\ to 
make the right side of eq. ([2]) equal zero. At the pole and the equator (9 = and 7r/2) we 
use the symmetric boundary condition: 

dpi 



89 
dQi 

~dT 

dv r 

~W 
v e = 0, 
ds i 

~d9 



0. 
= 0, 

0. 

0. 



(39) 
(40) 
(41) 
(42) 
(43) 



Due to the low Mach number of the expected flows, a direct c ompressible simulation 
is problematic, so adopting the same technique as iRempell (j2005bl ). we reduce the speed 



of sound by multiplying the right side of eq. ([T]) by 1/C 2 - The equation of continuity is 
therefore replaced with 



t£ + ^ div ^ v ) = °- 



(44) 



The speed of sound then becomes ( times smaller than the original speed. We use ( = 200 
in all our calculations. This technique can be used safely in our present study since we 
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only discuss stationary states, so t he factor ( beco mes unimportant. The validity of this 
technique is carefull y discuss e d by 



results presented by 



Rempell ( J2005bl ). We test our code by reproducing the 



Rempell ( I2005bl ) and check the numerical convergence by runs with 



different grid spacings. After checking and cleaning up at every time step, conservation 
of total mass, total angular momentum and total energy are maintained through the 
simulation runs. 



3. Stellar Differential Rotation and the Taylor-Proudman Theorem 



In this section, based on the work of 



Rempell (j2005bl ). we explain how the 



subadiabatically stratified region can generate solar-like differential rotation. The 
component of the vorticity equation can be expressed as 



du d 



[...] + r sin^ 



on 2 



g dsj 



(45) 



dt dz 7r 86 

where Q = Qq + ^i> and the z axis represents the rotational axis. The inertial term and 

the diffusion term are neglected. If the last term of eq. (1451) is zero, meaning there is no 
variation in entropy in the latitudinal direction, then dQ 2 /dz = in a stationary state, 
which is the Taylor-Proudman state. Solar-like differential rotation is generated in four 
stages. 

1. In the northern hemisphere, the A effect transports angular momentum in the negative 
z direction and generates a negative dQ 2 /dz. 

2. The negative dfl 2 /dz generates a negative u^ due to Coriolis force. This counter- 
clockwise meridional flow corresponds to a negative v r (downflow) at high latitudes 
and a positive v r (upflow) at low latitudes. 



3. As we mentioned in Section 12.21 downflow (upflow) generates positive (negative) 
entropy perturbations in the subadiabatically stratified layer beneath the convection 
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zone (5 < 0). Meridional flow can generate positive entropy perturbations at high 
latitudes and negative entropy perturbations at low latitudes. Therefore, ds\/d6 
becomes negative in the overshoot region. 

4. The negative ds\/d9 also keeps dfl 2 /dz negative in a stationary state. 

The profile of angular velocity in the convection zone is determined by a balance of 
angular momentum transport from meridional flow and a reduction in meridional flow from 
buoyancy force at the subadiabatic layer. 

4. RESULTS AND DISCUSSION 

We run simulations for seventeen cases, with Table [1] showing the parameters for each 
case. 

4.1. Stellar Differential Rotation 

In this section, we discuss the cases with angular velocities up to 16 times the solar 
value (represented by Qq), placing an emphasis on the morphology of stellar differential 
rotation. Fig. |3] shows the results of our calculations which correspond to cases 1-5 in Table 
[TJ It is found that the larger stellar angular velocity is, the more likely it is for differential 
rotation to be in the Taylor-Proudman state, in which the contour lines of the angular 
velocity are parallel to the rotational axis. To evaluate these results quantitatively, we 
define a parameter which denotes the morphology of differential rotation. We call it the 
Non- Taylor-Proudman parameter (hereafter the NTP parameter), which is expressed as 

"--mJ^-mK "^-^)^ (46) 
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where Qq is the angular velocity of the radiative zone. When the NTP parameter is zero, 
differential rotation is in the Taylor-Proudman state. Conversely, differential rotation is 
far from the Taylor-Proudman state with a large absolute value of the NTP parameter. 
The value of the NTP parameter with various stellar angular velocities is shown in Fig. HJ 
The NTP monotonically decreases with increases in stellar angular velocity. These results 
indicate that with large stellar angular velocity values, differential rotation approaches 
the Taylor-Proudman state. These results are counter-intuitive, however, since we do not 
expect differential rotation to approach the Taylor-Proudman state with increasing stellar 
angular velocity values, since the A effect, which is a driver of the deviation from the 
Taylor-Proudman state, is proportional to stellar angular velocity Qq. These are the most 
significant findings of this paper, so hereafter in this section we discuss these unexpected 
results. 

We next discuss the temperature difference between the equator and the pole at the 
base of the convection zone (r = 0.71i? Q ). Since temperature is given as a function of 
entropy by 

To 



7 



Si + ( T -1)- 
Po 



(47) 



and it is easier to measure than entropy, we use it here for discussing the thermal structure 
of the simulation results in the convection zone. Further, although it is mentioned 
in §3] that entropy gradient is crucial for breaking the Taylor-Proudman constraint, 
the temperature difference can be used as its proxy. Fig. \5\ shows the relationship 
between stellar angular velocity Qq and temperature difference AT at r = 0.71-R©, where 
AT = max(T 1 (rbc, 0)) — min(T 1 (rb c , &))■ Although the temperature difference monotonously 
increases with larger stellar angular velocity values, it is not enough to make the rotational 
profile largely deviate from the Taylor-Proudman state. This can be explained by using the 
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thermal wind equation, which is a steady state solution of eq. (j4"5j) : 



. n dQ 2 g dsi , An . 

= rsm Q -JL (48) 

oz 7r ou 

The inertial term and the diffusion term are neglected here. This equation indicates that, 
for a given value of the NTP, we need an entropy gradient proportional to Qq. However, 
our simulation results show that AT oc fif]' 58 , which means that as Qq increases, the 
thermal driving force becomes insufficient to push differential rotation away from the 
Taylor-Proudman state. In other words, the latitudinal entropy gradient in rapidly rotating 
stars is so small that differential rotation stays close to the Taylor-Proudman state. In our 
model, meridional flow generates latitudinal entropy gradient at the base of the convection 
zone. It is conjectured that the insufficient thermal drive is due to a slow meridional flow. 

We next investigate the dependence of meridional flow on stellar angular velocity. Fig. 
[6] shows the radial profile of latitudinal velocity vg at 9 = 45°, using the results of cases 1, 2 
and 9. In case 2, stellar angular velocity is twice that of case 1 (the solar value). In case 
9, stellar angular velocity is equal to the solar value, and the amplitude of the A effect is 
two times the value in case 1. Fig. [6] shows that meridional flow does not depend on stellar 
angular velocity, while it correlates with the A effect. Considering eq. (1M|) . the A effect 
increases with larger values of stellar angular velocity, since the amplitude of the A effect is 
proportional to Qq. The reason why differential rotation in rapidly rotation stars is close to 
the Taylor-Proudman state is that meridional flow does not become fast with large stellar 
angular velocity values. 

We interpret the result that the speed of meridional flow does not depend on stellar 
angular velocity in our model as follows. With large values of stellar angular velocity, more 
angular momentum is transported by the A-effect (Note that the A-effect is proportional to 
Qo in equation ( I34p ). so meridional flow obtains more energy from differential rotation. The 
energy gain does not result in an increase in speed because of the associated enhancement 
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of the Coriolis force, which bends the meridional flow in the longitudinal direction. Another 
explanation is possible in terms of angular momentum transport. The angular momentum 
fluxes from both meridional flow and the Reynolds stress (A effect) must be balanced in 
a steady state. The former is proportional to v m Qo and the latter is proportional to flo, 
where v m is the amplitude of meridiona 



on stellar angular velocity (IMiesch 



flow. Therefore, meridional flow does not depend 
20051 ). Our results (Fig. E]) indicate that with a larger 
stellar angular velocity (case 2), the above mechanism does not generate fast meridional 
flow. However, this does not occur when only the A effect is large (case 9). 

4.2. Angular Velocity Difference on the Surface 

In this subsection we discuss angular velocity difference AQ at the surface and 
the relationship between our results and previous observations. We conduct numerical 
simulations to investigate the physical process which determines AQ (cases 1, 6-11). We 
define angular velocity difference as AQ = max(f2i(r max , 9)) — min(fii(r max , 9)). 

AQ is determined by two opposing effects, a smoothing effect from turbulent viscosity 
and a steepening effect from the A effect. In a stationary state these two effects cancel 
each other out. Latitudinal flux for turbulent viscosity and the A effect can be written as 
PqVq w AVL/ A9 and po^oi^-o^O; respectively. Because these two have approximately the same 
value, AQ can be estimated as 

AQ ~ —A Q A9, (49) 

where A9 denotes the differential rotation region. 

In order to confirm eq. ( T4"9"j) . we conduct two sets of simulations, firstly varying the 
value of turbulent viscosity (fo v )> and secondly the amplitude of the A effect (A ). Note 
that the setting for turbulent viscosity does not reflect a real situation, since the coefficients 
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of turbulent viscosity and the A effect should have a common value. Nonetheless, this is 
necessary for the purpose of our investigation. The simulation results are shown in Figures 
[7] and [3 We obtain AQ oc ^°' 88 and AQ oc Aq 1 , which are consistent with eq. ( I49p . 

Fig. [9] shows the results of the dependence of AQ on Q (Cases 1-5). Asterisks denote 
the difference at the surface between the equator and the pole, squares show the difference 
between the equator and the colatitude 9 = 45°, and triangles are the difference between 
the equator and the colatitude 9 = 60°. The difference at low latitudes (squares and 
triangles) monotonically increases with stellar angular velocity. However this is not the case 
for angular velocity difference between the equator and the pole (asterisk). As we discussed 
in §4.1[ when stellar rotation velocity is large, the Taylor- Proudman state is achieved, 
meaning the gradient of angular velocity at the surface concentrates in lower latitudes. Due 
to this concentration, A9 becomes smaller in Eq. (I49J1 with larger values of Qq. Thus, AQq 
does not show an explicit dependence on Q . At low latitudes, A9 is fixed and the angular 
velocity difference increases with stellar angular velocity. We obtain AQ oc Q® A3 (between 
the equator and the colatitude 9 = 45°: squares) and Af2 oc Q®' 55 (between the equator and 
the colatitude 9 = 60°: triangles). This indicates that Af2/f2 decreases w ith stellar angu lar 



veloc i ty. These results are consistent with previo us stellar observations (IDonahue et al. 



1996 



Reiners &: Schmitt 



2003 



Barnes et al. 



20051 ) 



4.3. Variation of A-effect and superadiabaticity 

In this section, we discuss the dependence of meridional flow and differential rotation on 
free parameters. The parameter set is shown in Table [1] (cases 12-17). At first we investigate 
the influence of the variation of the A effect. The A effect has two free parameters, i.e., 
amplitude A and inclination angle A (see §2.4p . Amplitude is thought to become smaller 
with a larger stellar angular velocity, due to the saturation of the correlations such as (v' r v'^ 
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and (v'qv'A, where v' r , v' e and v'± are the radial, latitudinal and longitudinal component 
turbulent velocities, respectively. Fig. [H] shows that meridional flow becomes slower with 
a smaller A , keeping the f2 value constant (Case 10). It is clear with the result of 
§4.11 that meridion al flow becomes slow with a larger angular velocity when the variation 



of A is included. 



Brown et al 



( 120081 ) reported this effect with their three-dimensional 
hydrodynamic calculation. When meridional flow is slow, the entropy gradient generated by 
the subadiabatic layer is small, and differential rotation approaches the Taylor-Proudman 

state. 

The inclination angle is thought to be small with lar ge stellar angular velocity y alues, 



19931). In 



since the motion across the rotational axis is restricted (jKichatinov fc Riidigerl 
case 12, differential rotation with a small inclination angle (A = 2.5°) is calculated. Other 
parameters are the same as case 1. The radial distribution of meridional flow is shown in 
Fig. [61 Meridional flow becomes faster with a smaller inclination angle. Because of the 
efficient angular momentum transport in the z direction when the inclination angle is small, 
the second term on the right hand side of Eq. (T4"5l) is large. This generates a large u^, i.e. 
fast meridional flow. 

In summary, we found that rapid stellar rotation causes two opposing effects on the 
speed of meridional flow. The speed is reduced by the suppression of Ao, while it is enhanced 
by the angular momentum transport along the axial direction with a smaller A. Although 
the results of the three-dimensional calculation suggest that meridional flow becomes slower 
with a larger stellar angular velocity, our model cannot draw a conclusion about the speed 
of meridional flow in rapidly rotating stars. 

Next we investigate the influence of superadiabaticity in the convection zone. In 
cases 13-17, superadiabaticity in the convection zone 5 C = 1 x 10~ 6 . The differences of 
the NTP parameters with adiabatic and superadiabatic convection zones (-P n t P (<5 c =o) — 
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P ntp (,5 c= io-6))/-Piitp(5 c =o) are shown in Fig. [TUl The NTP parameter values with a 
superadiabatic convection zone are smaller than those with an adiabatic convection zone, 
since meridional flow in the super adiabatic convec tion zone makes the entropy gradient 
small. This result is suggested by iRempell ( 12005ai ). Note that the difference between the 



values of the NTP parameters with an adiabatic and those with a superadiabatic convection 
zone decreases as the stellar angular velocity increases, since the generation of entropy 
gradient by the subadiabtic layer becomes ineffective with a larger stellar angular velocity. 

5. SUMMARY 

We have investigated differential rotation in rapidly rotating stars using a mean field 
model. This work is significant because it can be used as a base for further research on 
stellar activity cycles, which are most likely caused by the dynamo action of differential 
rotation in the stellar convection zone. 

First, we investigated the morphology of differential rotation in rapidly rotating 
stars. Although more angular momentum is transported by convection with larger stellar 
angular velocity, the Coriolis force is stronger than in the solar case, so meridional flow 
does not be fast. In our model, meridional flow generates latitudinal entropy gradient in 
the subadiabatically stratified overshoot region. Since the meridional flow is not fast, the 
entropy gradient is insufficient to move differential rotation far from the Taylor-Proudman 
state in rapidly rotating stars. As a result, the differential rotation of stars with large stellar 
angular velocity is close to the Taylor-Proudman state. 

The temperature difference between latitudes is probably controlled by two important 
factors, i.e., the subadiabatic layer below the convection zone and anisotropic heat transport 
caused by turbulence and rotation. We suggest that the former is important in slow rotators 
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like the sun, and the latter in rapid rotators. The subadiabatic-layer effect is included 
in our model, while anisotropic heat transport is not. We found that the effect of the 
subadiabatic layer can generate a temperature difference AT = 10 K in the solar case, which 
moderately increases with higher rotation speeds, an d AT = 30 K in case Qq = 8Q Q . The 



three-dimensional simulations by 



Brown et al. 



(120081 ) include a self-consistent calculation 



of anisotropy of turbulent thermal transport but not the subadiabatic layer at the bottom 
boundary. In their calculation AT is most likely smaller than 10 K in the solar case, 
since they cannot reproduce the solar differential rotation only with anisotropy of thermal 
transport. Also, AT = 100 K in case f2 = 5f2 , which is larger than the case with the 
subadiabatic layer. We speculate that anisotropic heat transport becomes more significant 
in rapidly rotating stars. There is also a possibility that our calculated entropy gradient at 
the base of the convection zone can be used as a boundary conditio n for a self-consistent 



three dimensional simulation of stell ar convection (Miesch et al. 



rotation in rapidly rotating stars in 



20061 ). Note that differential 



Kiiker fc Stixl (120011 ) is not in the Taylor-Proudman 



state when anisotropy of turbulent thermal conductivity is included. A future study of 
the simultaneous effects of the attached subadiabatic layer beneath convection zone and 
anisotoropy of the turbulent thermal conductivity on stellar differential rotation would 
provide a better understanding of stellar differential rotation. 

Next, we investigated angular velocity difference at the surface. The A effect causes 
spatial difference in the rotation profile, while turbulent viscosity reduces the difference. 
Angular velocity difference AQ is determined in eq. ( H9|) . which is then used to investigate 
differential rotation in rapidly rotating stars. Since stellar rotation is close to the 
Taylor-Proudman state, and the radiative core is rotating rigidly, differential rotation is 
concentrated at low latitudes with large stellar angular velocity. This concentration leads 
to a small A8 in eq. ( )49l) . Therefore, only at low latitudes our model is consistent with 
stellar observations. 
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Our conclusions are as follows: (1) Differential rotation approaches the Taylor- 
Proudman state when stellar rotation is faster than solar rotation. (2) Entropy gradient 
generated by the attached subadiabatic layer beneath the convection zone becomes relatively 
small with a large stellar angular velocity. (3) Turbulent viscosity and turbulent angular 
momentum transport determine the spatial difference of angular velocity Afl. (4) The 
results of our mean field model can explain observations of stellar differential rotation. 

Our future work will focus on the stellar MHD dynamo. Several investi gations have 



been conducted on the stel 



200f 



Charbonneau fe Saar 



ar dynamo using a kinematic 



200f 



Moss fe Sokoloff 



2009 



dynamo framewo rk (JDikpati et al. 



Jouyeetal 



20101 ). Since, under 



such a framework, only the magnetic induction equation is solved using a given velocity 
field, solving a linear equation, such analysis does not give sufficient information on the 
strength of the dynamo-generated stellar magnetic field. To obtain the full amplitude of the 
stellar magnetic field, the feedba ck to the veloci ty field is required, i.e., an MHD framework. 



Adopting a similar approach to 



Rempell ( 120061 ). we can use the results of this paper to 



investigate the strength of the stellar magnetic field. Recent observations of the strength 
of the magnetic fi eld generated by stellar differential rotation have been conducted using 



spectroscopy (e.g. 



Petit et al. 



20081 ). A comparison of these observations and numerical 
calculations of the stellar dynamo could give new insight into the stellar magnetic field. 
Finally, our stellar MHD dynamo study would also contribute to t he understanding 



of recent investigations int o stellar magnetic cyclic activity periods ( INoyes et al 



Saar fc Brandenburg 



199a ). 
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Fig. 1. — Profiles of density, pressure and temperature as a function of radial distance in the 
reference state. This stratification is adiabatic. 
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Fig. 2. — Profiles of diffusivity as a function of radial distance for cases 1-5 and 9-11. The 
solid line denotes the coefficient of turbulent conductivity K t . The dashed line denotes the 
coefficient of the A effect v t \. The dash and three dots line denotes the coefficient of turbulent 

viscosity u tv . 
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Fig. 3. — Rotation profiles of the simulation results. Panels (a)-(e) correspond to cases 1-5, 
respectively. The stellar rotation rate for each case is given at the top of each panel. The area 
of red and solid lines (blue and dashed lines) rotates faster (slower) than the rigidly rotating 
core at the bottom boundary Color bars are given for angular velocity Q/27T = (Q + Qi)/2ir 
in the unit of nHz. The dotted lines in each panel indicate the base of the convection zone 
(r = 0.71i2©) and the colatitudes = 30° and 6 = 60°. 



-30- 



Q. 





NTP paramter 




i i i i i 




s 


10" 2 


r Xo i 








s 

N 


10~ 3 


r N o i 




s 




X ! 




s 




^ 
X 


10- 4 


r V i 








s 


10~ 5 


~i i i i i~ 



430 



860 1720 3440 6880 

^ /27i [nHz] 



Fig. 4. — NTP parameter as a function of stellar angular velocity VLq/2-k. The dashed line 
is the fit to the results showing a power-law function with an index of —2.9. 
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Fig. 5. — Temperature difference at the base of the convection zone (r = 0.71i? Q ) as a 
function of stellar angular velocity (Qq/27t). The dashed line is the fit to the results showing 
a power-law function with an index of 0.58. 
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Fig. 6. — Profiles of latitudinal velocity (vg) at colatitude 9 = 45° as a function of radial 
distance. In case 1, stellar angular velocity is the solar value, and the amplitude of angular 
momentum transport Ao = 1. In case 2, stellar angular velocity Qq — 2f2 Q . In case 9 and 
10, amplitude of the turbulent angular momentum transport Ao = 2 and 0.5, respectively. 
In case 12, the inclination angle of A effect A = 2.5°, and other parameters are the same as 
case 1. 
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Fig. 7. — Angular velocity difference at the surface as a function of the coefficient of turbulent 
viscosity z/ 0v . The dashed line is the fit to the results showing a power-law function with an 
index of -0.88. 
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Fig. 8. — Angular velocity difference at the surface as a function of the amplitude of the 
angular momentum transport A . The dashed line is the fit to the results showing a power- 
law function with an index of 1.1. 
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Fig. 9. — Angular velocity difference in three regions. Asterisks, squares and triangles 
represent the equator and the pole, the equator and the colatitude 6 = 45° and the equator 
and the colatitude 9 = 30°, respectively. The dashed and dotted lines are the fits to the 
results showing a power-law function with indices of 0.43 (squares) and 0.55 (triangles). 
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Fig. 10.— The difference of the NTP paramters between cases with adiabatic and supera- 
diabatic convection zone (P n tp(<$ c =o) ~ -Pnt P (5 c =io- 6 ))/-Pnt P (5 c =o) as a function of stellar angular 
velocity f2 /27r. The dashed line is the fit to the results showing a power-law function with 
an index of —1.4. 



